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ABSTRACT 


At  far-regional  and  near-teleseismic  distances  the  early  body-wave  coda  contains  information  that  is  potentially 
useful  to  monitoring  seismologists.  However,  waveforms  from  this  distance  range  are  typically  under-utilized 
because  of  propagation  complexities  that  cause  significant  difficulties  in  seismogram  interpretation.  For  example, 
the  first  ~20  seconds  of  a  far-regional  seismogram  often  include  multi-pathed  arrivals  caused  by  the  interaction  of 
the  wavefield  with  upper  mantle  discontinuities  at  220  km,  410  km  and  660  km  depth.  Depth  phases  (c.g.^pP, 
pP410,  sP,  sP4lO,  etc.)  also  add  complexity  to  the  early  part  of  the  seismogram  since  they  can  constructively  or 
destructively  interfere  with  the  primary  phase  arrivals.  Array  observations  from  earthquakes  in  central  Asia  regularly 
exhibit  a  variety  of  complex  phase  phenomena,  such  as  back-azimuth  anomalies,  emergent  or  late-arriving  first 
arrivals,  large  amplitude  secondary  arrivals,  and  interference  phenomena  between  upper  mantle  arrivals  and  depth 
phases. 

We  have  developed  array-based  methods  to  improve  the  characterization  of  primary  and  early  coda  phase  arrivals 
observed  at  far-regional  and  near-teleseismic  distances.  These  techniques  include  improved  signal  processing  to 
accurately  measure  the  delay  times  (r’s)  and  slownesses  (p’s)  of  primary  and  secondary  phases  from  small-aperture 
arrays.  We  use  these  r-p  measurements  to  develop  representative  crust  and  mantle  velocity-depth  profiles  and  suites 
of  synthetic  seismograms  through  those  models.  Then  we  use  the  processed  array  beams  to  derive  'wavefield 
templates';  i.e.,  grouped  observations  with  similar  phase  characteristics.  We  analyze  these  wavefield  templates  by 
comparing  them  with  synthetics,  looking  for  quantitative  explanations  for  the  phase  behaviors  we  observe. 

Our  approach  results  in  a  methodology  that  improves  phase  characterization  and  yields  earth  models  that  more 
accurately  predict  the  succession  of  expected  arrivals.  We  have  applied  our  techniques  to  observations  from  two 
regional  small -aperture  arrays  in  Kazakhstan  (MKAR  and  KKAR),  and  find  that  our  results  provide  insight  into 
wavefield  behaviors  that  are  regularly  observed  on  the  MKAR  and  KKAR  arrays. 
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OBJECTIVE 

The  objective  of  our  research  is  to  enhance  the  usefulness  of  regional,  small-aperture  arrays  by  developing 
array-based  methods  that  can  more  accurately  characterize  far-regional  (14"-29“)  seismic  wavefield  structure.  Far- 
regional  (14*-29“)  seismograms  are  particularly  interesting  because  of  the  strong  influence  upper-mantle 
heterogeneity  has  on  the  early  P  coda  arrival  structure  (<  20  secs).  The  suite  of  observed  phases  often  includes 
multi-pathed  arrivals  and  their  depth  phase  counterparts,  generated  from  wavefield  interaction  with  upper  mantle 
discontinuities  at  220  km,  410  km,  and  a  660-km  depth.  The  P-coda  arrivals  can  be  closely  spaced  in  time  (1-2  sec) 
and  have  similar  ray  parameters  that  range  between  8-1 1  sec/degree.  Further  wavefield  complexity  arises  from 
interference  between  depth-phase  multiples  and  near-receiver  scattered  arrivals  with  the  primary  arrivals.  These 
complexities  can  be  region  and  earthquake  specific. 

The  regional  seismic  arrays  that  have  been  built  in  the  last  fifteen  years  should  be  a  rich  data  source  for  the  study  of 
far-regional  phase  behavior.  The  arrays  are  composed  of  high-quality  borehole  seismometers  that  make  high- 
fidelity,  low-noise  recordings.  However,  beyond  regional  distances,  the  small  apertures  of  these  arrays  (<  5km)  limit 
their  usefulness  beyond  first-arrival  P-  and  S-  onset  picks.  At  these  distances  standard  array  methods 
(e.g.,  slant-stacking  and  frequency-wavenumber  analysis)  cannot  resolve  the  arrival  times  and  slownesses  of  primary 
and  secondary  arrivals,  making  the  proper  identification  and  classification  of  far-regional  arrivals  quite  difficult.  We 
are  overcoming  this  limitation  by  applying  refined  array  processing  techniques  in  conjunction  with  straightforward 
wavefield  generalization  methodologies.  Our  goal  is  to  characterize  the  commonly  observed  early  coda  arrivals  that 
propagate  from  the  different  seismic  regions  of  central  Asia,  utilizing  recordings  from  the  Makanchi  (MKAR)  and 
Karatau  (KKAR)  small -aperture  arrays  in  Kazakhstan  (Figure  1).  This  research  will  improve  the  usefulness  of 
small-aperture  arrays  by  increasing  their  ability  to  classify  small-magnitude  events  that  may  be  poorly  recorded 
regionally  and  teleseismically. 

RESEARCH  ACCOMPLISHED 

We  have  developed  array-based  methodologies  to  characterize  the  P  coda  of  far-regional  events  using  small  aperture 
arrays.  Our  methodology  includes  several  components:  1)  improved  small-aperture  array  processing  to  measure 
delay-times  (t)  and  slownesses  (p)  of  primary  and  secondary  arrivals;  2)  construction  of  region-specific  velocity- 
depth  profiles  derived  from  the  x-p  measurements  (i.e.,  x-p  transformations),  as  well  as  models  derived  from 
wavefield  continuation  methods,  and  3)  the  derivation  of ‘wavefield  templates’  from  clustering  of  array  beams  to 
capture  the  commonly  observed  arrival  structure.  These  methods  employ  both  theoretical  and  empirical  techniques 
to  reduce  a  large  waveform  data  set  (-600  earthquakes;  Figure  1)  into  a  smaller  set  of  robustly  observed  wave 
phenomena.  We  then  attempt  to  explain  these  phenomena  using  well-accepted,  straightforward  techniques.  The 
separate  components  of  our  methodologies  are  discussed  below. 


Figure  1.  Map  of  central  Asia  showing  the  location  of  the  Makanchi  (MKAR)  and  Karatau  (KKAR)  arrays 
in  Kazahkstan,  as  well  as  earthquakes  (circles)  used  in  this  study.  Colored  bands  mark  the  14’-29’ 
distance  range  from  each  array. 
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Improved  Small-Aperture  Array  Processing 

Our  objective  for  the  array  processing  is  to  compile  a  set  a  phase  delay  times  (t),  slowness  measurements  {p)  and 
array  beams  for  our  waveform  data  set  of  600  central  Asian  earthquakes.  The  small  aperture  (<  5  km)  of  many  of  the 
newer  arrays  installed  for  nuclear  monitoring  purposes  presents  a  serious  challenge  for  typical  array  processing 
methods.  However,  we  have  developed  and  applied  improved  velocity  spectra  analysis  (vespa)  and  phase-weighted 
stacking  methods  (PWS;  Schimmel  and  Paulssen,  1997)  to  improve  phase  identification  at  small-aperture  regional 
arrays.  Our  results  indicate  that  in  many  cases  we  can  detect  closely-spaced  arrivals  by  their  slowness  values  in  a 
time  window  that  includes  the  direct  P  arrival,  depth  phases  and  arrivals  from  upper-mantle  discontinuities 
(Ferris  and  Reiter,  2007). 

In  Figure  2  we  show  an  example  of  applying  phase-weighted  stacking  to  a  typical  regional  event  recorded  by  the 
MKAR  array.  At  a  distance  of  14.33“  and  depth  of  17  km,  a  series  of  arrivals  are  predicted  to  occur  within  the  first 
16  seconds  of  this  seismogram.  In  addition  to  the  main  P  branch,  the  reference  model  iasp91  (Kennett  and  Engdahl, 
1991)  predicts  the  pP  and  sP  depth  phases  at  a  delay  time  of  ~4  seconds  and  ~7  seconds,  respectively,  relative  to  the 
first  arrival.  These  depth  phases  are  followed  by  their  multi-pathed  counterparts  that  bottom  at  the  410  km 
discontinuity  and  have  a  delay  time  of  ~12  and  ~14  seconds  for  pP4I0  and  sP410.  When  compared  to  the 
decomposed  wavefield  (represented  by  the  phase  coherence  vespagram  in  Figure  2),  it  is  not  clear  how  these 
predicted  arrivals  correspond  to  the  observations.  Even  the  main  P410  arrival,  which  should  be  near  the  caustic  of 
the  travel-time  branch  and  thus  amplified,  is  not  obviously  identifiable.  However,  there  may  be  interference  between 
the  P4I0  and  sP  phase,  as  they  are  predicted  to  arrive  within  0.3  seconds  of  each  other.  If  the  earthquake  depth  is 
decreased  to  10  km,  the  pP  and  sP  arrivals  agree  better  with  the  observations.  However,  a  change  in  focal  depth  does 
not  account  for  all  the  observed  arrivals,  especially  the  coherent  later  arrivals  near  1 5  seconds,  or  explain  the 
ambiguity  in  the  P410  arrival. 


orld  yearday  time  Lat  Lon  z  distance.  1 4.33* 

000120;  2005243  124537.2  3637*  69.2r  17  km  seaz;  227.6’ 

—  -  *  —  -  ^  -  —  *  —  -  >  ■  ■  ‘  •  ■  ■  ■  -  ^  .  1  .  ■  ■  I  •  ■ 
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Figure  2.  Example  of  phase-weighted  array  processing  for  an  earthquake  recorded  by  the  MKAR  array  at  an 
epicentral  distance  of  14.3.  The  top  panel  shows  the  array  beam,  computed  using  the  slowness  of  the 
first  arrival.  Dashed  lines  mark  the  computed  onset  times.  The  middle  panel  shows  the  array  gather 
and  move-out  curves  for  detected  arrivals.  The  bottom  panel  shows  the  phase-weighted  coherence 
vespagram.  Black  circles  mark  arrivals  with  phase  coherence  greater  than  85%  (i.e.,  detections); 
black  squares  mark  the  computed  onset  time.  Colored  squares  show  the  iasp91  predicted  arrivals 
(yellow  =  P  branch  ,  orange  depth  phases,  green  =  sP  depth  phases). 
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This  exercise  of  comparing  our  array  measurements  to  those  predicted  from  global  reference  models  illustrates  the 
need  for  more  region-specific  models.  The  next  component  of  our  improved  phase  characterization  scheme  is  to  use 
our  array  measurement  and  array  beams  to  derive  models  that  can  more  accurately  match  the  observations.  Our 
methodologies  to  accomplish  this  and  preliminary  results  are  discussed  in  the  following  sections. 

Regional  Velocity  Depth  Profiles 

To  address  the  need  for  more  applicable  earth  models  and  to  generalize  our  array  observations,  we  are  developing 
two  separate  methods  to  derive  velocity-depth  profiles  for  the  specific  regions  monitored  by  the  MKAR  and  KKAR 
arrays.  Both  methods  are  based  on  x-p  (i.e.,  delay-time/slowness)  techniques  that  decompose  wavefields  into  their 
components  parts.  The  first  method  involves  directly  inverting  array-based  delay-time  and  slowness  measurements 
to  solve  for  a  velocity-depth  function  (e.g.,  Carbonell  and  Smithson,  1994).  The  second  method  applies  wavefield 
continuation  techniques  to  array-beam  record  sections  to  accomplish  the  same  goal  (e.g.,  McMechan,  1984).  While 
these  methods  are  not  new,  they  are  not  typically  applied  to  small-aperture  array  data,  but  rather  to  data  from  long 
linear  arrays  of  seismometers  (e.g.,  Morozov  et  al.,  2005).  Since  we  apply  these  methods  to  multiple  events  records, 
a  main  challenge  is  accounting  for  the  effects  of  differing  source  parameters  (origin  time,  event  depth,  and  focal 
mechanism),  which  may  be  unknown  or  poorly  constrained.  While  our  methods  are  under  active  development, 
preliminary  results  show  that  we  are  able  to  extract  reasonable  models  for  specific  regions  in  our  study  area.  The 
two  methods,  direct  inversion  and  wavefield  continuation,  are  discussed  below. 

Direct  Inversion  of  Array-Based  Delay-Time  and  Slowness  Measurements 

Figure  3  illustrates  how  we  prepare  array-based  delay-time  and  slowness  measurements  for  velocity-depth  inversion. 
In  this  example  we  groom  the  outliers  from  data  set  of  measurements  made  at  MKAR,  consisting  of  1 10  earthquakes 
ranging  from  14°-29°  epicentral  distance,  which  includes  earthquakes  extending  from  the  Hindu  Kush  region  of 
Pakistan/Afghanistan  to  the  Makran  coast  and  Zagros  Mountains  of  Iran.  The  grooming  removes  everything  outside 
the  8.0-14.5  (sec/deg)  slowness  range,  to  help  eliminate  measurements  from  P-io-S  scattered  signals  and  coherent 
noise.  We  also  discard  measurements  that  exhibit  significant  slowness  smearing  (i.e.,  measurement  uncertainties) 
during  the  PWS  processing.  Once  the  groomed  measurements  (Figure  3a)  have  been  collected,  we  reduce  them  to  a 
single  x-p  curve  by  computing  the  error-weighted-mean  r  value  within  evenly  spaced  slowness  bins  (Figure  3b), 
excluding  bins  with  only  a  few  measurements.  Averaging  the  r’s  within  a  common  slowness  bin  acts  to  correct  for 
errors  in  focal  depth  by  accounting  for  the  ±x  offset  between  P  and  pP.  While  this  is  not  a  perfect  correction  due  to 
missed  phases  in  the  array  processing,  the  inclusion  of  sP  arrivals,  and  measurement  error,  we  have  found  it  more 
effective  and  feasible  than  correcting  x  on  an  event  basis  using  catalog  depths,  which  typically  have  large 
uncertainties.  The  averaged  x-p  curves  are  the  observations  that  are  input  into  the  velocity-depth  inversions.  The 
curve  shown  in  Figure  3b  generally  agrees  with  the  iasp91  x-p  curve;  however,  it  is  not  as  smooth,  which  may  be 
partially  caused  by  measurement  error  and  inadequate  averaging,  not  just  differing  earth  structure. 

Figure  4  shows  a  preliminary  velocity  profile  determined  from  the  MKAR  x-p  data  shown  in  Figure  3.  There  are 
many  ways  to  parameterize  this  inverse  problem,  and  the  results  in  Figure  4  represent  one  solution.  To  derive  the 
model  we  solved  for  the  thicknesses  of  a  sequence  of  layers  that  best  fit  the  x-p  observations,  using  a  damped  least- 
squares  formulation.  For  this  particular  inversion,  we  linearized  the  problem  by  assigning  each  slowness  p  to  the 
layer  representing  its  maximum  depth  of  penetration  (i.e.,  the  ray  turning  point),  which  depends  on  the  initial 
velocity-depth  profile  and  is  updated  after  each  inversion  iteration.  The  under-determined  system  of  equations  we 
invert  is  r  =  G  /i,  where  r  is  a  vector  of  our  zero-offset  time  observations  and  h  is  vector  of  model  layer 
thicknesses.  Each  row  of  G  is  composed  of  the  vertical  slowness  through  each  layer  up  to  the  turning  depth  for  each 
p  observation.  While  linearizing  the  inversion  allows  a  direct  estimation  of  model  uncertainties,  the  final  results  tend 
to  be  strongly  influenced  by  the  initial  model.  We  are  currently  examining  other  algorithms  to  better  understand 
model  non-uniqueness,  as  well  as  other  parameterizations  of  the  linear  inverse  problem. 
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Figure  3.  Sample  r-p  data  from  the  MKAR  array,  a)  The  groomed  MKAR  data  (blue  circles);  red  dots  show 
the  iasp91  r-p  curve  for  a  surface-focus  source.  The  range  in  r  for  a  particular  slowness  is  caused 
by  the  difference  in  r  for  earthquakes  at  different  depths,  rather  than  just  measurement  error,  b) 
The  averaged  T-p  curve  from  the  groomed  data,  evenly  sampled  in  slowness  at  0.15  sec/deg. 
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Figure  4.  Preliminary  T-p  inversion  for  a  1-D  velocity  profile  of  the  MKAR  data  shown  in  Figure  3.  a)  final 

model  (red)  compared  to  the  starting  model  (black);  b)  the  model  resolution  matrix;  and  c)  the  data 
resolution  matrix. 


Wavefield  Continuation  Analysis  using  Array  Beams 

Wavefield  continuation  (e.g.,  Clayton  and  McMechan,  1981)  is  commonly  used  to  extract  velocity  profiles  from 
record  sections  recorded  by  long  linear  arrays  of  sensors.  In  our  case,  array  beams  form  the  record  sections.  We 
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arrange  the  beams  according  to  their  epicentral  distance  and  align  them  based  on  the  first  arrivals  (Figure  5a). 
Arrivals  are  aligned  relative  to  a  reference  model,  a  time  shift  that  acts  as  an  event-based  static  correction.  The 
correction  tries  to  account  for  the  differences  in  earthquake  depth  and  origin  time.  This  is  a  simplistic  correction  and 
does  not  correctly  account  for  the  time  shifts  of  secondary  phases  (i.e.,  pre-critical  arrivals  from  the  major 
discontinuities),  whose  arrival  times  relative  to  the  first  arrival  depend  on  earthquake  depth.  However,  some  type  of 
correction  needs  to  be  made  to  each  array  beam  in  order  to  get  a  reasonable  x-p  transformation.  To  account  for  the 
difference  in  source  mechanism,  we  perform  the  wavefield  continuation  calculations  on  the  envelopes  of  the  beams 
rather  than  the  raw  amplitude  records.  This  acts  as  a  low-pass  filter,  reducing  the  resolution  of  the  resulting 
velocity-depth  profile. 

The  wavefield  continuation  method  consists  of  two  linear  transformations  of  the  wavefield  record  sections.  First, 
record  sections  are  slant  stacked,  transforming  them  from  the  distance-time  domain  into  the  x-p  domain  (Figure  5b). 
This  is  followed  by  a  downward  continuation,  or  depth  migration,  of  the  wavefield  to  transform  the  x-p  data  to  the 
slowness-depth  plane.  This  second  step  is  an  iterative  process,  where  the  x-p  data  are  repeatedly  migrated  until  the 
slowness-depth  image  converges  to  the  input  velocity  model  (Figure  5c).  There  are  several  ways  to  update  the 
velocity  model  between  iterations.  We  are  currently  using  a  scheme  that  picks  the  maximum  amplitude  at  each  depth 
from  the  slowness-depth  image  and  then  computes  the  weighted  average  between  it  in  the  input  velocity  model.  This 
scheme  performs  adequately;  however,  noise  in  the  initial  x-p  transformation  can  cause  artifacts  in  the  final  velocity- 
depth  profile.  We  are  currently  working  on  methods  to  improve  the  iterative  migration  procedure. 
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Figure  5.  Example  showing  the  transformation  of  a  record  section  (a)  to  the  x-p  domain  (b)  and  the 

downward  continuation  to  the  slowness-depth  domain  (c).  Array  beams  are  from  104  earthquakes 
that  extend  from  the  Hindu  Kush  region  of  Pakistan/Afghanistan  to  the  Makran  coast  and  Zagros 
mountains  of  Iran. 
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Figure  6  shows  a  comparison  of  our  preliminary  velocity-depth  profile  derived  from  the  wavefield  continuation 
procedure  to  the  iasp91  model.  Below  400  km  depth,  our  model  exhibits  moderately  slower  velocities  than  iasp91, 
and  a  gradient  zone  at  the  410  km  rather  than  a  sharp  discontinuity  (Figure  6a).  A  comparison  of  the  slowness  vs. 
distance  curves  for  the  predicted  P  arrivals  for  each  model  are  shown  in  Figure  6b.  Here  the  differences  are  more 
apparent,  where  our  derived  model  shows  more  structure  above  the  410  km  discontinuity  (slowness  <  ~12  sec/deg), 
and  a  significant  difference  in  the  caustic  points  of  both  the  410  km  (termination  of  B-C  branch  in  Figure  6b)  and 
660  km  discontinuity  (termination  of  D-E  branch).  While  our  model  is  preliminary,  some  of  the  features  we  observe 
may  be  related  to  an  inadequate  distance  sampling  of  our  array  beams,  as  well  as  the  windowing  procedure  used  in 
the  initial  x-p  transformation.  Both  of  these  would  explain  the  initial  appearance  of  the  predicted  410  km  and  660  km 
triplicate  arrivals  at  greater  distances  (18°  and  20°,  respectively)  in  our  derived  model  than  the  iasp91  model,  where 
they  appear  at  14°  and  18°,  respectively. 
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a)  b) 


Figure  6.  Comparison  of  the  velocity-depth  profile  from  wavefield  continuation  to  the  iasp91  reference 
model:  (a)  derived  model  (black)  and  iasp91  model  (red),  (b)  Slowness  vs.  distance  comparison 
between  the  derived  model  (black)  and  the  iasp91  model  (red)  for  predicted  P  arrivals. 

Waveform  Templating  for  Groups  of  Like  Waveforms 

Another  component  of  our  phase  characterization  methodology  is  the  analysis  of  the  array  beams  using  clustering 
techniques.  In  addition  to  the  arrival  times  and  slowness  of  the  coda  phase  arrivals,  phase-weighted  stacking  analysis 
yields  array  beams,  which  are  then  input  to  a  waveform  clustering  algorithm.  The  algorithm  is  based  on  ‘fuzzy 
clustering’  to  form  groups  of  beams  that  exhibit  similar  characteristics;  it  has  been  used  widely  in  pattern 
recognition  applications  (Bezdek,  1981).  Clusters  are  defined  by  assigning  each  beam  a  cluster  membership  value, 
which  is  a  quantitative  measure  that  incorporates  a  distance  measure  between  each  array  beam  and  a  representation 
of  the  cluster  centers.  We  have  experimented  with  distance  measures  involving  LI  and  L2  norms,  as  well  as 
waveform  semblance,  with  varying  degrees  of  success  that  depend  on  data  noise  levels.  From  each  cluster  of  similar 
beams,  we  derive  a  template  waveform  using  a  stacking  process  that  weights  each  cluster  member  by  its  degree  of 
membership  to  that  particular  cluster.  The  ‘fuzziness’  of  the  method  allows  a  single  beam  to  be  a  member  of  more 
than  one  cluster  group.  The  objective  behind  the  waveform  clustering  is  to  reduce  the  database  of  observations  to  a 
set  of  representative  waveforms  that  exhibit  consistent  phase-arrival  behavior.  Sometimes  the  clustering  is 
geographic,  but  there  are  other  wave  phenomena  that  can  also  produce  groupings  of  similar  waveforms 

Figure  7a  shows  the  cluster  group  and  its  associated  wavefield  template  that  results  from  applying  our  waveform 
clustering  algorithm  to  the  set  of  100  events  shown  in  Figure  7b.  These  beams  are  derived  from  seismograms 
recorded  by  the  MKAR  array  in  the  14°-29°  distance  range  that  extends  from  northern  Pakistan  to  the  Makran  coast 
in  Iran.  In  this  example,  the  algorithm  clusters  the  events  into  5  groups,  with  12-24  members  per  group.  We  applied 
an  L2-norm  distance  measure  and  aligned  the  array  beams  prior  to  clustering  using  multi-channel  cross  correlation 
(VanDecar  and  Crosson,  1990),  applied  to  the  first  four  seconds  of  each  beam  signal.  We  have  found  the  alignment 
of  the  beams  to  be  beneficial  to  the  clustering  procedure.  However,  the  length  of  the  cross-correlation  window  is  an 
important  variable.  For  example,  if  the  window  is  too  long,  the  array  beams  align  on  the  maximum  amplitude  signal 
in  the  record.  For  some  beams,  the  maximum  amplitude  arrival  may  occur  early  in  the  seismogram,  while  for  others 
it  occurs  much  later.  This  can  result  in  some  unusual  clustering  of  events  that  seems  counter-intuitive.  For  this 
reason,  we  find  clustering  methods  based  strictly  on  cross  correlation,  as  used  in  other  types  of  studies  (e.g.,  Menke 
1999;  Ferretti  et  al.,  2005;  Barani  et  al.,  2007)  are  not  suitable  for  our  particular  needs. 
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Figure  7.  Example  of  wavefield  clustering  to  generate  a  wavefield  template,  (a)  The  bottom  panel  shows 
the  cluster  members  sorted  by  epicentral  distance  (shown  on  the  left).  The  green  band 
highlights  the  time  band  of  cross  correlation.  The  top  panel  shows  the  computed  wavefield 
template  (black)  and  l-o*  deviation  (red),  (b)  A  record  section  of  all  100  earthquakes  input  to  the 
clustering  algorithm. 


For  the  example  shown  in  Figure  7,  the  resulting  cluster  members  include  earthquakes  that  are  near  the  14.3° 
distance  range,  but  some  more  distant  earthquakes  are  also  included.  The  earthquake  depth  range  is  also  variable, 
but  may  be  due  to  poorly  constrained  catalog  depths.  In  general,  the  cluster  members  show  consistent  waveform 
structure  in  the  cross-correlated  portion  of  the  signal.  Later  arrivals  show  more  variability  between  the  cluster 
members,  but  still  appear  to  align  well  at  approximately  8  seconds  and  later.  The  most  variability  occurs  between 
4  and  8  seconds  delay  time,  where  some  records  show  large  amplitude  arrivals  while  other  show  little  signal.  Some 
of  this  variability  is  likely  due  to  differences  in  earthquake  distance,  triplicated  arrivals,  and  the  presence  of  depth 
phases.  We  continue  to  refine  and  improve  the  wavefield  clustering  algorithm. 

Phase  Characterization  Analysis 


To  characterize  observed  phase  arrivals,  we  calculate  different  misfit  measures  between  our  suite  of  body-wave 
synthetic  seismograms  and  the  observed  wavefield  templates.  Since  the  synthetic  waveforms  are  generated  from 
models  derived  from  the  r-p  data,  arrival  times  of  the  wavefield  templates  should  be  well  matched.  We  note  that  to 
derive  meaningful  measurements,  we  must  account  for  frequency  content  and  amplitude  variations. 

Figure  8  provides  an  example  of  how  we  fit  the  synthetic  seismograms  to  the  wavefield  template  shown  at  the  right 
of  Figure  7.  We  computed  synthetics  using  a  test  model  (CRUST2.0  over  AK135)  for  focal  depths  ranging  from 
0-60  km  and  event  distances  between  14°  and  30°.  Then  we  calculated  the  L2  root-mean-square  (RMS)  residual  of 
the  differences  between  the  wavefield  template  and  each  synthetic  seismogram.  The  left  panel  in  Figure  8  shows  the 
results  of  this  exercise,  with  a  black  cross  denoting  the  distance  (15.4°)  and  depth  (32  km)  at  which  the  minimum 
RMS  value  occurs.  In  general,  the  best-fitting  synthetic  fits  the  overall  structure  of  the  template,  even  though  the 
method  is  fairly  crude.  However,  a  portion  of  the  wavefield  template  does  not  fit  well  (between  7-8  s),  and  this  may 
be  due  to  a  number  of  causes.  For  instance,  the  test  model  we  used  to  generate  the  synthetics  may  be  inaccurate. 
Also,  the  focal  mechanism  or  the  misfit  measure  we  employed  in  the  comparisons  could  be  inadequate. 
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Figure  8.  Example  of  fitting  synthetic  seismograms  to  an  empirical  wavefield  template.  Right  panel  shows  the 
wavefield  template  (black)  and  the  best-fitting  synthetic  (red).  The  left  panel  shows  an  image  of  the 
root  mean-square  residual  between  the  wavefield  template  and  a  suite  of  synthetics  computed  for 
varying  event  distance  and  depth.  The  best-fitting  synthetic  for  this  wavefield  template  corresponds 
to  a  distance  of  15.4^  and  an  event  depth  of  32  km  (black  cross). 


CONCLUSIONS  AND  RECOMMENDATIONS 

Regional  seismic  arrays  that  have  been  recently  installed  for  nuclear  monitoring  are  under-utilized  in  the  study  of 
far-regional  arrivals.  The  small  aperture  of  many  of  these  arrays  (<  5km)  restricts  their  usefulness  at  these  distances 
beyond  first  arrival  onset  picks  of  P  and  S  waves.  However,  our  research  is  overcoming  this  limitation  by  applying 
refined  array  processing  techniques  in  conjunction  with  x-p  and  wavefield  templating  analysis  methods.  Our 
methodology  improves  the  characterization  of  primary  and  early  coda-phase  arrivals  observed  at  far-regional  and 
near-teleseismic  distances.  Our  approach  is  to  distill  the  wide  variety  of  seismicity  we  observe  to  subsets  of 
commonly  and  robustly  observed  phase-arrival  phenomena.  We  are  then  using  well-accepted  modeling  and 
inversion  techniques  to  explain  these  phenomena.  Our  aim  is  to  explain  as  much  of  the  phenomena  as  we  can  with 
simple  and  straightforward  techniques,  leaving  the  anomalies  for  future  research. 

We  are  developing  and  applying  our  techniques  to  central  Asian  earthquakes  recorded  on  the  MKAR  and  RKAR 
arrays  in  Kazahkstan.  Our  results  indicate  that  we  can  differentiate  between  the  numerous  arrivals  of  the  early  coda. 
However,  global  reference  models  cannot  capture  the  phase  succession  and  arrival-time  behavior  we  observed  from 
the  complex  tectonic  regions  of  central  Asia.  To  address  this,  we  are  developing  regional  1-D  models  directly  from 
the  array  measurements  (delay-times,  slowness,  and  beams).  Since  these  models  are  derived  from  the  data,  they 
should  better  explain  the  phase  behavior  we  observe  from  specific  regions.  To  test  these  models  and  further  improve 
phase  characterization,  we  are  constructing  ‘wavefield  templates’  through  cluster  analysis  to  generalize  the 
waveform  structure  from  the  different  seismic  regions.  The  templates  are  then  used  in  a  waveform  fitting  analysis  to 
gain  a  better  understanding  of  the  phase  phenomenon  observed  from  the  complex  seismic  regions  of  central  Asia. 
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